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Abstract: We present an analysis of Bose-Einstein condensation for a system of non- 
interacting spin-0 particles in a harmonic oscillator confining potential trap. We 
discuss why a confined system of particles differs both qualitatively and quanti- 
tatively from an identical system which is not confined. One crucial difference 
is that a confined system is not characterized by a critical temperature in the 
same way as an unconfined system such as the free boson gas. We present the 
results of both a numerical and analytic analysis of the problem of Bose-Einstein 
condensation in a general anisotropic harmonic oscillator confining potential. 
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1 Introduction. 

One of the most interesting properties of a system of bosons is that under certain 
conditions it is possible to have a phase transition at a critical value of the temperature 
in which all of the bosons can condense into the ground state. It is now well over seventy 
years since the phenomenon referred to as Bose-Einstein condensation (BEG) was first 
predicted for the ideal nonrelativistic Bose gas 0. Nowadays it is well known to 
happen if the spatial dimension D > 3. (See for the case D = 3 and @] for general 
D.) 
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Until recently the best experimental evidence that BEC could occur in a real physi- 
cal system was liquid helium, as suggested originally by London [^]. However although 
the behaviour of liquid helium at low temperatures can be qualitatively described by 
the free boson gas model, the detailed behaviour deviates substantially from this sim- 
ple model. Physically this is of course because the effects of interactions which are 
neglected in the free boson gas model are important in liquid helium. More recently it 
was sugg ested i, that BEC could occur for excitons in certain types of non-metallic 
crystals (such as CuCl for example). There is now good evidence for this in a number 
of experiments p|. 

An important development in the last year has been the experimental attempts 
to observe BEC in very cold gases of rubidium , lithium [0 , and sodium |n| . 



This experimental work has stimulated theoretical studies to try to understand the 



underlying physics of the situation [12, 13, 14]. The systems are very dilute and as 



a first approximation would be expected to be well described by a boson gas model 
with no interactions among the atoms. The atoms are confined in complicated magnetic 
traps which can be modelled by harmonic oscillator potentials. There have been several 



studies of BEC in harmonic oscillator confining potentials ||T5|, |T^, |T^, 0, |T^. The 
purpose of our paper is to examine the condensation of bosons in a harmonic oscillator 
potential in a detailed way which does not use the density of states approach of Refs. [|1^, 
[17| , [l^, |l9l. Unlike the situation for a boson gas with no external confining potential 



in free space there does not exist a critical temperature which signals a phase 
transition. However, we will show that there is a temperature at which the specific 
heat has a maximum which can be identified as the temperature at which BEC occurs. 
(A short report of our results was given in Ref. [pl[] .) 

The fact that a gas of bosons (neglecting interactions) in a harmonic oscillator 
potential does not have a phase transition at some critical temperature is already 
apparent from the early work of [|l5l. (This will also be shown below in a very simple 
way.) A similar situation occurs for a system of charged bosons in a homogeneous 
magnetic field; in three spatial dimensions BEC does not occur in the same way as for 
the case where there is no magnetic field E^. The same is true for bosons confined 



by spatial boundaries [^. (See also Ref. p4|.) A general criterion to decide whether 
or not BEC occurs has been given recently by us ||2^, and it is easy to show that the 
criterion is not met for a system of bosons confined by a harmonic oscillator potential. 
Since the experiments of Refs. |10|, ^ are claiming to observe BEC, a natural 



question which arises concerns the exact nature of the phenomenon. If BEC as found 
normally for the free boson gas is impossible for a system of bosons in a confining 
potential, then in what sense does BEC occur ? A natural criterion which has been 
used in systems of finite size has been to look at the maximum of the specific heat . 
We will apply this criterion to the case of bosons confined by a harmonic oscillator 
potential. Given the details of the harmonic oscillator potential trap, it is possible to 
calculate a characteristic temperature which can be compared with the values found 
in the experiments. 

The paper is organized as follows. In section 2 we consider a gas of bosons in an 
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isotropic harmonic oscillator potential. Although the potentials present in the exper- 
iment are anisotropic, we start with this model because it is much easier technically 
and as we will see afterwards, the anisotropy has not much influence on the critical 
temperature where the specific heat has its maximum. Thermodynamical quantities 
are given in terms of sums resulting from the grand partition function of the system. In 
section 3 we present the technique of how to obtain approximate analytic results for all 
sums involved. These techniques are used in section 4 to derive simple expressions for 
the thermodynamical quantities which allow for a determination of the critical temper- 
ature and the groundstate occupation number. Combining our analytical techniques 
with numerical calculations, we present the detailed behaviour of the specific heat and 
other quantities of interest. In section 5 we generalize our model to the example of an 
anisotropic oscillator. Once more the detailed behaviour of several quantities relevant 
for the recent experiments is given. Appendices A-D contain several technical details 
on how to obtain the approximate results for all quantities of interest. Finally, the con- 
clusions summarize all important results, give a brief comparison between our method 
and that of Refs. 0, 113 , and present a short discussion of further work. 



2 The isotropic harmonic oscillator potential 

For mathematical simplicity let us start with the case of an isotropic harmonic oscillator 
potential. We will assume that the system can be described by a grand canonical 
ensemble. The grand potential is defined by 

q = -Y,ln{l-zexp{-(3E^)) , (2.1) 

N 

where /3 = {kT)~^, En are the energy levels, and z = e^'^ is the fugacity in terms of the 



chemical potential fi. It proves convenient to expand the logarithm in (|2.1|) to obtain 



oo 



q=Y.-Y.^M-nl3EM). (2.2) 

n=l ^ N 

For an isotropic harmonic oscillator characterized by an angular frequency uj the energy 
levels are given by En^n2nz = hjj{ni + + + 3/2), ni,n2,n3 e Nq. If we set 
^1 + ^2 + ^3 = k, where k G Nq, then the energy levels may be ordered in the way 
Ek = {k + ?,/2)hu with multiphcity {k + l){k + 2)/2. The sum over N in (U) may be 
performed to obtain 

oo n/3(fj,-3/2huj) 

g=E ^n_.-n.^3 > (2-3) 



n=l 



n(l - e-"^)3 ' 



cles is given by N = P ^ ( ) , which becomes 



where we have defined the dimensionless variable x = huj/{kT). The number of parti- 
ng ^ 

(2.4) 



n=l 



3 



when ( |2.3| ) is used. 

In order that the number of particles remains positive, it is necessary for fi < 
{3/2)hu. (More generally, we require fi < Eq where Eq is the lowest energy level.) 
Normally the critical temperature for BEC is the temperature at which n = Eq, which 
for the isotropic harmonic oscillator reads fi = {3/2)hu!. It is now easy to see that BEC 
cannot occur in the same way for bosons confined in the harmonic oscillator potential 
as it does for bosons in free space. In the case of the free boson gas in free space with 
no confining potential, as the temperature is lowered the chemical potential fi increases 
from negative values towards the value 0. (This is in agreement with the general result 
fi = Eq quoted above since the lowest energy level is zero for the free boson gas.) The 
value of the temperature at which fi = defines a critical temperature Tc determined 
in terms of the particle density. At temperatures lower than T^, fi remains frozen at 
the value fi = 0, and the number of particles found in excited states is bounded. If the 
total number of particles exceeds this bound then the only possibility is for the excess 
particles to be found in the ground state, giving rise to BEC. This standard scenario 
is described in |^D| in some detail. The phase transition which occurs is related to 
the breaking of the U{1) gauge symmetry associated with the change of phase of the 



Schrodinger field. We have discussed this in a recent review P7 |. 

The origin of the different behaviour between the confined and the free Bose gas 
might be seen more clearly by considering the number of particles in the ground state 
with energy {3/2)huj. In addition to the dimensionless quantity x we introduce fi = 
huj{3/2—e), the limit e — corresponding to the limit of the chemical potential reaching 
its critical value. In terms of x and e, the number of particles in the groundstate is 

-Aground T ■ (^•^) 

gex — i 

For e — > we have N ground —>■ oo, this being essentially the reason that no BEC in the 
usual sense that e reaches at some finite temperature might occur. For given x and 
particle number it is clear from (|2.5|) that e > [l/x) ln((A^ + 1)/N) and that e can 
reach only in the zero temperature limit or in the limit —>■ oo. However, as is also 
clear from (p.5|) , for fixed particle number A^, once e is small enough, it is essentially 
only the groundstate which is occupied. Lowering the temperature, which is the same as 
increasing x, this will happen unavoidably as can be seen from ( |2.4| ). The temperature 
at which the groundstate starts to increase considerably its occupation number is the 
most dramatic moment in the system. The extreme behaviour is similar but not equal 
to a phase transition. Although quantities are changing rapidly, everything behaves 
smoothly; no discontinuities appear. The argument presented here may be performed 
in a much more general context using the setting described in [^. Further discussion 



of this important point we postpone until we have presented our analytical analysis of 
the thermodynamical quantities. 

Let us continue with the internal energy U of the system which is given in terms of 
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the grand potential q by 



U 



d /i (9 



In terms of the dimensionless quantities x and e, U reads 

U_ _ _dq _ (3/2 - e) dq 
huj dx X de 

Using the series for q given in ( |2.3|) results in 

U 3 



hu 2 



where 



N + 3ui 



:i-e 



-nx\—4: 



(2.6) 



(2.7) 



(2.^ 



(2.9) 



n=l 



In terms of U, the specific heat reads 



AT.t^ held fixed 



(2.10) 



Since N is held fixed when computing C, only ui contributes to the specific heat, and 
we find 



C/k = -3x 



2 I ^"^1 



dx 



N,ui held fixed 



(2.11) 



Once more due to fixed N, alternatively one might use 



oo ^—enx 



Ui 



E 



(2.12) 



„^,(l-e— )'' 

which differs from Ui only by a multiple of A^. Continuing with (|2.11|) one first finds 



dui 
dx 



-So 



N,uj held fixed 



1 + g^iexW 



-4^3, 



with 



^ne-"^™(l-e" 



-nx\—A 



and 



n=l 



oo 
n=l 



:i-e- 



(2.13) 



:2.14) 



:2.i5) 
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Differentiating equation ( |2.4| ) with respect to x for fixed N,uj, — — is determined. 

— ox 

We find 

V ox /jv,.. held fixed 
where in addition to 5*2 and S3 we introduced 

= ^ne-"^"(l-e-"")"^- (2-17) 

n=l 

(2.18) 

So putting the results of equations ( p.l3| ) and ( |2.16|) together, we arrive at 



C/k = |4^^ gSj . (2.19) 

Before we proceed with the numerical analysis of some of the above thermody- 
namical quantities, we turn now to the analytical treatment of these quantities for 
some range of parameters x and e. Let us look at the relevant range of parameters in 
the sodium experiment. (Qualitatively it will be the same for the other experiments.) 
There |Tl[] we have u/{27r) = 416 Hz if we use the geometric mean of the frequencies. 
The relevant temperature range is around 2fiK. It is therefore seen, that the behaviour 
of the thermodynamical quantities for small x is desired. A plausible approach is to 
argue that for x <^ 1, it is justified to replace sums which have arisen in the expansions 
above with integrals. Care must be exercised with this to take a proper account of the 
density of states. (We will return to this in Sec. 6.) This is tantamount to regarding 
the energy levels as continuous rather than discrete. However, we have shown recently 
that the behaviour of thermodynamical systems with a discrete energy spectrum is 
completely different from one with a continuous energy spectrum. In the first case, 
no real BEC can occur whereas in the second case it does (By real BEC, we 



mean that there is a phase transition such as that which occurs in the free boson gas.) 
If the correct behaviour for small x is desired, one approach which is definitely safe 
is to deal with the exact sums. The sums do not converge very rapidly for small x, 
nor do they display in any transparent way the behaviour at small x. However, it 
is possible to convert the sums into contour integrals, and by deforming the contours 
of integration in an appropriate way obtain at least asymptotic expansions for some 
appropriate range of the parameters. The details of this procedure will be described 
in the following section. 
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3 Analytical treatment of harmonic oscillator sums 



Let us now describe in some detail the analytical treatment of the sums appearing in 
the thermodynamical quantities. As one can see, all sums involved are of the form 

oo —nex—mnx 

f{l,k,m) = T r, (3.1) 

with different integral values for l,k,m. Let us consider the range of x <^ 1, which 
is actually fulfilled in the recent experiments described below, and in addition e <^ 1 
corresponding to the range where a phase transition occurs in 3 dimensional free space 
and, as we will see, corresponding to the range very close to the maximum of the 
specific heat. 

Probably the best technique for the analysis of (|3.1|) in the mentioned range of 
parameters is the use of the Mellin-Barnes integral representation. We are going to 
apply it in its simplest form, making use of 



^ c+joo 



2m 

c—ioo 



da r(a)w-", (3.2) 



valid for 3fJf > and c G R, c > 0. Equation ( |3.2| ) is easily proven by closing the 
contour to the left obtaining immediately the power series expansion of exp(— w). In 
order to apply equation ( ^.2[ ) write ( |3.1| ) in the form 



oo 



f{l,k,m) = Y^n-^ e-'^^[^+'^+^^+-+^*] 

n=l ii,...,ifc=0 

c+ioo 



oo oo 1 / 

= ^ c/ar(a)n-"x-"[e + m + ii + ...+ifc]~"(.3.3) 

n=l n,...,*,=o 

Now we would like to interchange the summation and integration in order to arrive at 
an expression in terms of known zeta functions, in detail the Riemann zeta function 



n=l 

and a Barnes zeta function ||28|| , 

oo 

CB{s,a,k)= Y [h + ■■■ + ik + a]~' . (3.5) 

«lvi«fc = 

The basic properties of the Barnes zeta function are summarized in the Appendix A. 
In order to allow for the interchange of summations and integration one has to ensure 
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the absolute convergence of the resulting sums |30|. This is certainly true for 
> max {k, 1 — /) and one arrives at 



c+ioo 



f{l,k,m) = — : / dar{a)x "(^{0 + l)(is{a,m + e, k). (3.6) 

c—ioo 

This is a very suitable starting point for the analysis of certain properties of the sums 
/(/, fc, m). Closing the contour to the right corresponds to the large-x expansion; closing 
it to the left to the small-x expansion. To the right of the contour the integrand in 
( |3.6|) has no poles, which means that the large-x behaviour contains no inverse power 
in X. One might show however, that the contribution from the contour itself is not 
vanishing at infinity leading to exponentially damped contributions for x ^ 00, the 
well known behaviour of partition sums at low temperature. 

As mentioned, this is not the range of interest for recent experiments and we con- 
centrate on the small-x behaviour thus closing the contour to the left. Closing to the 
left there appear to be three different sources of poles, 

i. ) poles of Ce(a, m + e. A;) for a = 1, k] (See Appendix A.) 

ii. ) pole of Cr{(^ + /) for a = 1 — /; 

iii. ) poles of r(a) for a = —p, p G Nq. 

Depending on the value of / there might be a double pole at a = 1 — /. In detail for 
I = 0, —1, 1 — k, there is a double pole from i.) and ii.); for / = 1, 00 there is a 
double pole from ii.) and iii.). Collecting all poles is an easy exercise for all values of 
/. We will restrict to the relevant values of / for our problem, these being I = 1,0, —1. 
For Z = 1 we find 

k 

f{l,k,m) = ^r(n)x-Xi?(l + ^)^es CB(^,e + m,A;) (3.7) 

n=l 

m + e,k)- (lna;)CB(0, e + m,k)+ 0{x), 

Res (b being the residue of the Barnes zeta function, and its derivative with respect 
to s. (See equation (^.5p.) The residues and values of the Barnes zeta function are 
derived in Appendix A. The leading important residues are 



Res Csik, a, k) 



1 



{k-l)V 



k — 2a 

Res CB{k-l,a,k) = ~ 2)V ^^'^ 

6a^ -6ak + {k/2)(3k-l) 



Res (sik — 2, a, k) 



12(A;-3)! 



Also the derivative of Res Cb('S? k) with respect to s at s = might be determined 
in terms of derivatives of the Hurwitz zeta function, but the contributions are of sub- 
leading order and will not be used here and thus are not presented. 
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For / = 0, —1, the formula analogous to ( p.7| ) reads 



k 

f{l,k,m) = ^ T{n)x~'^(R{l + n)Res (i3{n,m + e, k) 



n = l 



+x'-^r(l - /) {PP Ce(l -l,m + e,k) 

+ (7 - Inx + - l))Res (ei^ -l,m + e,k)} 
+CRil)CB{0,m + e,k) + O{x), (3.9) 

PP (is denoting the finite part of Cs, ip{x) = (d/dx) Inr(x) and ip^l) = —7. 

The presented asymptotics together with ( |A.7| ) allow one to obtain the small x and 
e behaviour of the theory. In the next section we list the asymptotic expansions of the 
various thermodynamical quantities. The needed sums are given in the Appendix B 
for the convenience of the reader. 

By changing slightly the procedure described previously, it is also possible to obtain 
an expansion for x <^ 1 not restricted to e ^ 1. To find this representation write instead 
of equation ( |3.3| ) 

00 -nxe 00 1 c+«oo 

f{l,k,m) = Y.^- TT- / ^^«r(«)n-"x-"[m + 2i + ...+z,]-", (3.10) 

n=l ii,...,ifc=0 

which is found by using equation (|3.2|) only for exp{—nx[m+ii + . . .+ik]) and excluding 
the part exp(— nxe) from the procedure. For the case that m = 0, one has to treat 
separately the zero-mode ii = . . .ik = 0. Closing the contour once more to the left to 
obtain the x <C 1 behaviour, an expansion in terms of the polylogarithm 

00 I 

L^n{x)=Y.-. (3.11) 
1=1 ' 



is found. The basic properties of the polylogarithm may be found in |3^. For 
reasons of clarity the corresponding results are summarized in Appendix D, however 
only for the anisotropic harmonic oscillator presented in section 5. The isotropic case 
is easily extracted from Appendix D. Because our intention is to try to present sim- 
ple analytical expressions we will not use the expansion in polylogarithms, since by 
necessity this would involve numerical evaluation. 



4 Temperature dependence of the thermodynami- 
cal quantities 

Having described in detail the application of Mellin-Barnes integral techniques for the 
approximate calculation of harmonic oscillator sums, we are now prepared to present 
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the results for all thermodynamical quantities of interest. First of all for the grand 
potential we have g = /(I, 3, 0) and using equation (|3.7|) results in 

Cfi(4) , (i - ^) CniS) ^ Ci?.(2) , ,,,, , . 
q = — 5 \- \ \- 0{\nx,\ne). (4.1) 

3j Olj Ob 

For the number of particles we have = /(O, 3, 0) and find 

^^c«M^ (i-)<"P) ,l^o(!^). ,4.2) 

x-^ ex \ X 



It is also possible to present equation ( |4.2|) in a slightly different way. As we have 
explained in some detail in section 2, a quantity of special interest is the number of 
particles in the groundstate, given by ( |2.5| ). Splitting = N ground + N^xcited and using 
the same techniques as described in section 3, but now with = . . . = ifc = excluded 
from the summation in equation (|3.3| ) (this summation index actually corresponds 
exactly to the groundstate), the following asymptotic expansion is found, 

,4,3) 

Equation ( ^4.2|) is very useful to determine e as a function of N and x. This then leads 
with the help of equation (|4.3|) to the groundstate occupation number as a function of 
X (for fixed N). 

Using ui = /(0,4, 1) and furthermore equation ( p.8|) , the asymptotic expansion of 
the internal energy reads 

U 3Cr(4) ^ US) /9 _ , 13Cij(2) 



huj x^ x^ \2 J Ax"^ 

The asymptotics for ui together with the asymptotics of 5*1, 5*2, 5*3 needed for the 
analysis of the specific heat are listed in appendix B. After some calculation we find 
using ( |2.19| ) the following result, 

C 12C/j(4) 9Ci;(3) 2Cn{2) 12eU3) 

+ —T- + — -3 (4-5) 



186^Ci^(2)Cfi(3) _ 96^Ci?(3)^ ^ 96^Ci^(2)Ci^(3)^ 

OC^ 



+ O {\nx, ^ 



This concludes the list of the asymptotic behaviour for x <^ 1, e <^ 1, of the most 
important thermodynamic quantities considered here. 

The corresponding expansions in terms of the polylogarithm (see the end of section 
3) are also easily obtained and given in Appendix D. Once more the results for the 
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isotropic harmonic oscillator follow immediately from those for the anisotropic oscilla- 
tor. 

The above results together with a numerical treatment of the thermodynamical 
quantities using their explicit representations in form of the sums given in section 2 
allows a very detailed prescription over the whole temperature range of relevance for 
the recent experiments. As an illustration we will first choose parameters pertinent 
to the rubidium experiment p. We choose = 2000 here. The aim is to compute 
the chemical potential which is given by e. (Recall that fi = huj{3/2 — e).) This may 
be done by solving ( p.4|) for e as a function of x. (Recall that x = Tioj/ikT) is the 
inverse temperature in appropriate dimensionless units.) The numerical result of this 
calculation is shown as the solid curve in Fig. 1. As can be seen from this figure, e 
undergoes a very rapid decrease from values of the order of unity to values of the order 
of 10~^ over a very small range of x. After this sharp decrease, e goes asymptotically 
to zero as x increases ( or as T decreases). This contrasts with a real phase transition 
such as occurs in the free Bose gas where e would reach the value e = at some nonzero 
temperature which could be identified with the critical temperature. From ( ^.5| ) it can 
be seen that this sudden drop in e is associated with a sudden rise in the ground state 
occupation number, and is therefore associated with the onset of BEC. 

Although the chemical potential has a sudden change, the change happens in a 
completely smooth way; thus the identification of a specific critical temperature is 
problematic in this case. One approach which has been used in finite volume systems 
where similar behaviour occurs p6| is to calculate the maximum of the specific heat and 



identify the temperature at which the maximum occurs with the critical temperature. 
In Fig. 2 we illustrate with a solid curve the result of a calculation of the specific heat 
using the exact harmonic oscillator sums. It is seen to have quite a sharp but smooth 
maximum at a value Xm — 0.0921. If we use uj/{2'k) = 60 Hz, as for the strong trap 
referred to in Ref. |jl2|, then the specific heat maximum occurs at the temperature 



T ~ 3.127 X 10-8 K. 

We now turn from a numerical evaluation to the use of our approximate analytical 
results detailed above. Because our approximation assumed that x and e were both 
small we would not expect the results to apply for x < x„i — 0.0921 since Fig. 1 shows 
that e is already becoming quite large. If we define / to be the fraction of particles in 
the ground state, 

Nground = fN , (4.6) 



then from (|2.5|) we have 



ex = \n(l + ^y (4.7) 



Using ( [4.3| ) this yields 



;i - f)N ^ Cr{^)x-' + (^ - 6)C«(2)x-2 . (4.^ 
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e may be eliminated from ( [4.8| ) by using ( |4.7| ). This results in a cubic equation which 
determines x for a given / and A^. Once x has been found e is determined by ( [4.7|) . 
We have shown the result of this approximate evaluation of e as the diamonds in 
Fig. 1. As expected, once x decreases below Xm — 0.0921, the agreement between our 
approximation and the exact result breaks down. However over the region where the 
specific heat maximum occurs, our approximation for e is quite good. As x increases 
(so that the temperature becomes less that the critical temperature), the agreement 
between our approximate value for e and the true value becomes better and better. 
Increasing values of x correspond to an increasing fraction of particles in the ground 
state. Fig. 1 shows that the agreement between our approximate value for e and the 
true value remains remarkably good even up to values of e ^ 0.5. Given that we 
assumed e <^ 1 in our derivations, this is an unexpected result. 

Since our approximate result for e is good in the region where the specific heat 
maximum occurs, we can have some faith in our other approximate results for x > 
Xm — 0.0921. In Fig. 2 the diamonds illustrate the result of using our approximate 
result ( [4.5| ) for the specific heat. Again the agreement between the approximation and 
the exact result is seen to be good up to the maximum. For values of x < Xm, the 
approximation breaks down for the reason already mentioned. We have shown a more 
detailed comparison between our approximation for the specific heat and the true value 
in Fig. 3. The diamonds illustrate the ratio of our result to the exact value. For x ~ Xm 
our approximate result is within a few percent of the true value, and for x > Xm the 
agreement becomes better than 1%. 

We can also compare our results to the bulk results obtained by directly converting 
the harmonic oscillator sums into integrals as in Refs. [0, |19[. (We will use the ter- 



minology bulk rather than thermodynamic limit as in Ref. P^.) For the specific heat 



this amounts to just keeping the first term on the right hand side of ( [4. 51 ) : 

C,uik/k = 12Cii(4)x-=' . (4.9) 

For the particle number the bulk result consists of dropping the term in in ( ^.3| ) 
along with all subdominant terms, taking 

Nbuik = N ground + CR(3)a;-^ . (4.10) 



A way of improving the bulk approximation was given in Refs. ||T7| , p!8| , and we will 
return to the relationship of this improvement to our approach in Sec. 6. 

The bulk transition temperature is obtained by equating Nground to zero. This gives 



•^bulk 



(4.11) 

where Xbuik = ft-^ / (kThunS) . Eq. ([4.9| ) holds for x > x^uik- We have plotted the ratio 
of Chuik to the exact specific heat as the crosses shown in Fig. 3. Although the agree- 
ment between the bulk value and the exact value is quite good near the specific heat 
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maximum, it is off by about 25% when x reaches the value 0.4. In contrast our ap- 
proximation is off by less than 1%. Another feature of using the bulk approximation is 
that the specific heat is found to be discontinuous at the bulk temperature, in contrast 
to the smooth behaviour found in Fig. 2. 

A final comparison we will make is for the ground state occupation number. In 
Fig. 4 the diamonds represent the result of using our approximation for the particle 
number and the crosses the result of using the bulk value. Both results are shown as a 
ratio with the exact value found from a numerical evaluation of the harmonic oscillator 
sum. Close to the specific heat maximum our results are off by about 10% and the 
bulk results are off by around 300%. As x increases, our result converges very rapidly 
towards the true value, whereas the convergence of the bulk results is slower. For large 
X both approximations become indistinguishable. 

It is possible to obtain an approximate analytic expression for the BEC temperature 
and compare it to the bulk temperature. Suppose that x is close to the bulk value given 
in (|4.11|) and write 

X = Xbuikil + V) , (4.12) 
for some small i]. If we assume e << 3/2, then from ( |4.8| ) we find 



V 



f) 



lf+ ^^(^^ iV-V3 
3^^2[Ck(3)]2/3 



(4.13) 



This result assumes that Xbuik is small, but makes no assumption about the size of /. 
Typically we find that for particle numbers N ~ 10^ — 10^, at the point where the 
specific heat maximum occurs / is a few percent. We can therefore say (1 — f)~^ ~ 1, 
which leads to 



T-T, 



bulk 



T 



lf+ iV-V3 
3^ ^ 2[C«(3)]2/3 



(4.14) 



If we put / = this gives the result in Ref. |T^, |T8|. For = 2000 it is easily seen 
that using the bulk value for the temperature is only about 10% higher than using the 
value determined by the maximum of the specific heat. Both temperatures approach 
one another as the particle number increases. Thus as an estimate of the critical 
temperature, the use of the bulk value is quite a good approximation. However as our 
results above show, care must be exercised in using the bulk values for the particle 
number or specific heat. 

To finish our discussion of the isotropic harmonic oscillator we wish to mention 
briefly some results for other choices of particle number. We would expect that as A^ 
increases not only will the bulk approximation become better, but so will ours. We 
have done the calculations just described for N = 2x 10^, 2 x 10^ and 5 x 10^. Because 
the resulting figures are similar to those already presented in the case A^ = 2000 there 
is little point in showing them here. The expectation of improved agreement between 
the approximations and the exact result is borne out. The specific heat maximum 
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occurs at x„r ^ 0.0408 for = 2 x 10^, at Xm ^ 0.0185 for = 2 x 10^ and at 
Xm — 0.0136 for = 5 X 10^. These values correspond very closely to those obtained 
using the approximate formula ([4.14 ). 



5 The anisotropic harmonic oscillator potential 

After having described in detail the isotropic harmonic oscillator potential let us explain 
the new technical problems one encounters when treating the anisotropic case. The 
calculation parallels very much the one for the isotropic harmonic oscillator and we 
can be brief. The energy eigenvalues are given by 

Enmzns = h'^Ui (rii + ^\ni e Nq. (5.1) 
i=i ^ 

As we will see in the following, it is useful to introduce the dimensionless quantities, 

1 3 

Xi = hjSuOi] Vt = - '^uji\ a = %(3Vt. 
3 i=i 



We then have 



3 3 

i=l 



In analogy to the harmonic oscillator we use furthermore 

'3 
,2 



= hVt ( - — e 



to find 



3 

i=l 



In terms of the dimensionless variables the grand potential reads 



For the particle number we have 



N=y^^ r- (5.4) 



With eqs. ( |5.3| ) and ( |5.4| ) one easily gets the internal energy. 

U 3 
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where we defined 



•J^i \ — > J. 
Vj' — — / 

' a n?=i (1 - e-^^") (1 - e-"^') ' 



(5.6) 



Continuing for the specific heat as done for the analysis of the isotropic harmonic 
oscillator, we arrive at 



with 



1 



2A 



i=l 



-E 



X 



5-1 



1=1 

3 3 

i=l 1=1 



(5.7) 



and 



^2,1 



E 



ne 



inti(i-e-"^") 



- E 



ne 



in 



S3,ii — Yl 



ne 



—n€a—n{xi+xi) 



{ U.j=i (1 - e-^^") (1 - e-"^<) (1 - e-=^'") 

The asymptotic expansions of all above quantities may be obtained using the same 
techniques as for the harmonic oscillator described in section 3. The only difference is 
that one has to deal with the shghtly more general function 



CB(s,a|f) = ^(a + mf)' 



(5.8) 



m=0 



The asymptotic expansions for the thermodynamical quantities involves the residues 
of the function (^g(s,a|r), the basic properties of which are summarized in Appendix 
A. Using once more the Mellin-Barnes integral representation in complete analogy 
to section 3, we arrived at the asymptotics for Ui, Si, 82,1 and Sz,ii- These are all 
summarized in Appendix C. We here list only the asymptotics of the physical quantities, 



Cr{A) , CR(3)a ^3 _ ^ 



X1X2X3 X1X2X3 \2 
(^(2)0;^ ^ 3:1X2 + X1X3 + X2X3 3 
12q;2 4 



(5.9) 



X1X2X3 



V 



+ ..., 
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N = i£PL + ^«(H)^f!-.)+^ + ..., (5.10) 



X1X2X3 \2 ) ta 
U 3Ch(4) ^ 3Ch(3) /3 ^ (j^^j, 



X1X2X3 \2 72 ) 
C 12C«(4) ^ 9Cij(3)a 9a2e2^^(3)2 



/C ^1X2X3 XiX2X^ {XiX2Xz] 



2 



4 X1X2X3 12 X1X2X3 
12e«Ci?(3) 18a2e2(^(2)C^(3) 



X1X2X3 (a;iX2X3)2 

4^4a /'o^/" ^'Q^2 



9aVCfi(2)Cfi( 3) 

(xiX2X3)3 



The above asymptotic expansions are found to be a good approximation close to, but 
lower in temperature, the maximum of the specific heat. Once more we are able to 
present a numerical as well as an analytical calculation of the relevant quantities. 

Suppose that we define / to be the fraction of particles in the ground state as we 
did in Sec. 4. Nground is given by 

ground 1^ 

From ( |5.10| ), noting that the term in l/(ea) arises from the ground state, we find the 
number of particles in excited states is given by 

^.,,.M!)5!±^0i(f!(^^^^^)(|_,)J^. (5,3) 

ujiUj2ijJz a'^ o ^co'iCi;2 "^i^^s ijJ2^?, V2 J 

We will illustrate the results for the case of = 2000 as for the isotropic harmonic 
oscillator using the frequencies for the rubidium experiment. (The results shown are 



independent of whether the strong or weak trap is used, since the difference is one 
of an overall scaling of the oscillator frequencies, and the results we use are independent 
of such a scaling.) Fig. 5 shows the result of a comparison of our approximate result 
for e (illustrated with diamonds) and the exact result illustrated by the solid curve. 
Again the agreement is quite good even for relatively large values of e. Fig. 6 shows 
the comparison between our approximation for the specific heat in ( p.l2| ) and the 
exact value found from the harmonic oscillator sums. The maximum occurs for a ~ 
0.106. Fig. 7 shows the ratio of our approximation to the exact specific heat, and 
for comparison the result of using the bulk expression. As for the isotropic oscillator 
calculations, the bulk expression shows a significant deviation from the true result; 
when a ~ 0.2 the bulk result is off by about 15%, whereas our approximation is within 
about 1% of the true value. Fig. 8 shows the ratio of the approximate particle numbers 
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to the true value. Our result is seen to have better agreement close to the specific heat 
maximum, but both our result and the bulk result rapidly converge towards the true 
value as a increases. 

We can now see that the anisotropy has only a small effect on the critical tem- 
perature. Using the frequencies ojx = uj^ = 2407r/A/8 s~^, and ojj, = 2407r s~^, with 
a ~ 0.106 we find ~ 3.09 x 10~® K as the temperature at which the specific heat 
maximum occurs. This can be compared with the temperature of 3.13 x 10~^ K found 
in the isotropic case in Sec. 4. The bulk temperature (eq. ( |4.11|) still holds in the 
anisotropic case with u the geometric mean of the frequencies) is about 3.41 x 10~^ K. 



6 Conclusions 

In conclusion, in this article we presented a detailed analysis of several thermodynam- 
ical quantities for a system of non-interacting spin-0 particles in a general harmonic 
oscillator confining potential trap. Although there is no phase transition such as that 
occurring in the free unconfined boson gas, it is possible to identify a temperature at 
which BEG occurs by looking at the maximum in the specific heat. We have seen that 
this temperature is nearly identical to the temperature where the groundstate occupa- 
tion starts to increase considerably, the effect actually seen in the recent experiments 
through the peak in the velocity distribution of the sample. 

Attempting to compare the results obtained here directly with the experiments 
must be done with a certain degree of caution. In the first place we have ignored 
interatomic interactions, so that there is no distinction made between gases with a 
positive scattering length |TT|, and those with a negative scattering length |T^. 
Secondly, it is perhaps not quite so clear that the use of the grand canonical ensemble 
is justified for systems with such a relatively small number of particles [^]. With these 



caveats in mind, for the case of rubidium we found the specific heat maximum to occur 
at T ~ 31 nK for 2000 particles. For the case of sodium, if we use N = 2.5 x 10^ we find 
the specific heat maximum to occur at T ~ 1.16 /iK. The results for the temperature 
found from using the bulk approximation were very close to these values, so it is 
unlikely that the present experiments can distinguish between the bulk approximation, 
our approximation or the exact value. It was apparent from our calculations that the 
specific heat found from our approximation was much closer to the exact result than the 
bulk approximation was. Perhaps in future experiments it will be possible to provide 
a more stringent test of the various approximations. 

We now wish to mention the comparison between our results and the density of 
states method used in Ref. [|l^, |T^. In this approach the authors separated off the 
ground state contribution for the particle number and treated the remaining terms in 
the sum by approximating it with an integral over the energy. The density of states 
was parametrized by 
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where uj = (ujiuj2UJ3Y^^ and 7 is a dimensionless function of the frequencies. In the 
isotropic case 7 = 3/2, but the authors [|l^, had to determine 7 numerically in the 



anisotropic case. The bulk approximation we mentioned earlier consists of ignoring the 
term in 7. By contrast, the approximate method we presented is entirely analytical 
with no numerical evaluations required. By comparing the results of our calculation 
for the particle number with those of Ref. |jl^, |18[ we can deduce an analytic value for 
7. We find 

7 = -a;M + + . 6.2 

2 \UJ1UJ2 ^2^3 oosUJiJ 



This value has also been borne out by an independent evaluation ||3^ of the density of 
states which in addition obtains the next order correction to ( pTl| ) in the anisotropic 
case. 

Although the gases in the experiments are dilute, in order to get a quantitatively 
more satisfactory picture for the recent experiments, the vapor has to be treated as a 
weakly interacting system. In the quantum field theory approach to BEC this might 
be done in a systematic way [0. One possibility is to treat the interaction as a 



perturbation and calculate the leading corrections to the free boson gas treated in the 
present article. The detailed knowledge of the lowest order in perturbation theory 
provided here is thus an important basis for future developments in this direction. 
Another possibility is to consider an effective theory for the groundstate, its occupation 
number being a very good indicator for the onset of BEC. 

After this paper was submitted for publication, another independent calculation of 
BEC in harmonic oscillator confining potentials appeared ^5|. This paper uses the 
Euler-Maclaurin summation formula to evaluate the harmonic oscillator sums. The 
result is expressed in terms of polylogarithms, much like the results found in our Ap- 
pendix D. One difference between this approach and ours is that the authors of Ref. 
introduce an effective fugacity which means that the argument of their polylogarithms 
never becomes equal to unity. It is straightforward to modify the approach we have dis- 
cussed in Appendix D and obtain in an easy way results which appear to be equivalent 
to those of Ref. . However this method necessarily involves a numerical evaluation of 
the polylogarithms, as we mentioned earlier in connection with our own polylogarithm 
results, and is therefore not entirely analytic. 
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A The Barnes zeta function 

As we have seen, in order to determine the asymptotic expansion of some thermody- 
namical quantities, we need several properties of the Barnes zeta function 



CB{s,a\r)= ^(a + mr)-^ (A.l) 

m=0 

with r a d-dimensional vector. In equation (|3.5|) we used the notation Cb{s, a, d) for 
r = 1. The residues of Cb('S, a\r) at s = 1, d and the values of the function at s = — p, 
p G No are most easily deduced using the representation as a contour integral 

' ' ^ 27r Jc ^ ' nti(l -e-"»*) 

where the contour C is counterclockwise enclosing the positive real axis. The only 
possible pole occurs at t = 0. For that reason one might like to introduce the generalized 
Bernoulli polynomials |^ through 



^ ^ i:B^\a\r-)^-^. (A.3) 



In terms of these it is immediate that for = 1, 



[n — l)!(fl 

and for p G N, 



Res Cs{n, a\r) = , B^^tW (A.4) 



a r 



For the leading asymptotics of the thermodynamical quantities we need at most the 
first three residues in ( [A.4|) . These are given explicitly by 

Res Ceid, a\r) 



Etir.-2a 

2{d-2y.utin 



ResCs{d-l,a\r) = ,,r=^J^, , , (A.6) 



Res (s{d — 2, a\r) 



12{d-3)\UUr^ 
In addition to the above equation we needed only 

CB{s,a\r) = ^ + Oia'), (A.7) 

which is obvious from the original sum ( [A.l|) . Using equations ( |A.6| ) and ( |A.7|) we 
found the asymptotic expansion for all thermodynamical quantities. 
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B Asymptotics for x ^ 1, e ^ 1 of the sums ui, Si, 

S2 and 5*3 

In this appendix we list the results used for the derivation of the internal energy and 
the specific heat of the isotropic harmonic oscillator confining potential. We needed 
the following asymptotic expansions : 

ui = /(0,4,1) 

Ck(4) , {l-e)U3) , 1Ck(2) , ^fe I 



= /(-1, 3,0) 

1 ^ a(2) ^ ^(\rix\ 



S2 = /(-l,4,i: 



0(3) ^Ck(2) ^^/lnx^^\ ^ 

^3 = /(-1, 5, 2) 

Ci?(4) Ci^(3) 6C«(3) 1 Ci?.(2) , ^flnx e\ 
x5 2x4 x4 12 x3 ' 

These results lead, after some calculation, to equations (WM and ([4. 51). 



C Asymptotics for x ^ 1, e ^ 1 of the sums Ui, Si, 

S2,i and Ss^tj 

In this appendix we give the results used for the derivation of the asymptotics of several 
thermodynamical quantities. (See equations (|5.9|) - (|5.12|) .) 

First of all we need the analogous results to equations ( |3.7| ) and ( |3.9| ) for the 
anisotropic oscillator. Unfortunately for the anisotropic oscillator it is not possible 
to write all needed sums in a unified form as done in equation (p.l|). For that reason 
we have to list several results. The techniques are exactly the same techniques as those 
employed in section 3. We found for / = —1, 

3 

T{m)(ji{m + l)Res Csij^i ea\x) 

m — 1 

+r(l-/) {PP CB(l-/,e«|x) 

+ (7 + ^(1 - /)) Res Cb(1 - /, ea\x)] 
+Ci?(OCB(0,e«|x) + ..., 



oo 



{n^ nj=i (1 - e-^^") 
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whereas for Z = 1 one has 



3 



n=l '''lli=l \^ ^ ■' ) m=l 

+CB(0,ea|f) + ... . 

These results can be used for the calculation of q, N and 5*1. 
For Ui and 5*2, i one needs 

oo ^—nta—nxi 



\ X{j=i (1 - e-^^") (1 - e-"^*) 

4 



m=l 



+r(l -0{PP CH(l-^,e« + a;i|(f,Xi)) (7 + ^(1-0) 
xi?es — /, ea + Xi\{x^ ^i))} + ■■■ ■ 

Finally for S^^ij we need 



E 



1 n/=i (1 - e-^'^O (1 - e-"^') (1 - e-"^j ) 
5 

^ r(m)Cii(m - l)i?es Ce(?7i, ea + + | (f , Xi, xj)) 



m — 1 

m/2 



+PP Cb(2, ea + + I (f , Xi, xj)) 
+Res Cb(2, ect + Xj + Xj\{x,Xi,Xj)) + ... . 

These results are enough to find the following expansions, 

xiX2X^a 2x1X2X3 \ a 

^ CR(2)a _ _^ a:ia:2 + xix^ + 3:2X3 + ' 
12x1X2X3 y q; 
_ 1 , Cii(2) , 

(aej^ X1X2X3 
= Cii(3) ^ Cii(2)« X, 

XiX2X3Xj 2XiX2X3Xj V a 



S3,ij 



Ci?(4) , Ci?(3)a /3 Xi + Xj \ CR(2)a 



X]^X2^3^i^j X\X2X^X{Xj \2 2q; / Y'Zx \X 2X '^^X iX j 

/ Xj + X,- X1X2 + X1X3 + X2X3 + X^ + X^ + SXjX,- \ 

X 9 - 9 H ^ - 

\ a / 
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D Asymptotics for x <ti I for the anisotropic har- 
monic oscillator 



As mentioned in section 3 and 4, it is possible to obtain an asymptotic expansion valid 
for X <^ 1 without restricting e to the range of small parameters. The way how to 
obtain the approximation is described at the end of section 3. The results analogous 
to equation (p.7|) and ( p. 91) read 

n=l lli=l l-^ ^ ) 1=1 



oo —l—nea—nxi 4 

^ — ^ it^ fc- ^ — A 



T{n)Res CB{n,Xi\{x,Xi))Lii+ri (e + ... , 



n=l 



ne 



-nat—n{xi+Xj) 



(D.l) 



- e-"^«)(l - e-"^0 nf=i (1 - e~"^0 
^ T{n) X iies C,B{n,Xi + a;j|(x, Xj, (e"'") + ... . 

n=l 

As a result, the following asymptotics for the thermodynamical quantities are derived, 



1 



X1X2X3 ^ ' 2X1X2X3 

/3 1 X1X2 + X1X3 +X2X3\ 

+ + — ]Li2(e )+... 

V 4x1X2X3 12 X1X2X3 



-^L.3(e-)+^^L.2(e- 



X1X2X3 ^ ^ 2X1X2X3 

^ /3 ^ 1 X1X2 + X1X3 + X2X3 \ / _^ 

\4xiX2X3 12 X1X2X3 / 

m2 0x1X2X3 ^ ^ 2x1X2X3 ^ ^ 

+ -7; — (36q;^ + X1X2 + X1X3 + X2X3)Lz2 (e"'°) + 

120x1X2X3 ^ ^ 

k X1X2X3 [ Li2 (e ^ 

« r / X 27 Lzi(e-)Lzl(e-) l 



X1X2X3 [ ^ ^ 2 Lii (e-^°) 

This concludes the list of asymptotic expansions which we are going to give in the 
present article. 
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Figure Captions 



Figure 1 This shows e as a function of x = fiuj /{kT). The sohd curve shows the result 
found from a numerical evaluation of the exact harmonic oscillator sums for the 
isoptropic case for N = 2000. The diamonds show the result of using our analytic 
approximation. 

Figure 2 The specific heat computed numerically for the isotropic harmonic oscillator 
is shown as the solid curve. The diamonds show the result using our approxi- 
mation. The units for the specific heat are in factors of the Boltzmann constant 
k. The particle number is N — 2000. The maximum occurs for x ~ 0.0921. 
Our approximation breaks down for x below the point where the specific heat 
maximum occurs. 

Figure 3 The diamonds show the ratio of our approximation for the specific heat to 
the exact result. The crosses denote the ratio of the bulk specific heat to the 
exact value. N = 2000 is taken. Both results become increasingly inaccurate 
below the specific heat maximum. 

Figure 4 The ratio of the approximate to the exact ground state particle number is 
shown. The diamonds illustrate the result of using our approximation and the 
crosses denote the results found using the bulk approximation. 

Figure 5 The numerical result for e is plotted against a = h{uJi + uj2 + uj^) / (3/cT) 
for the anisotropic oscillator for = 2000. The frequencies arc taken to be 
wi/(27r) = W2/(27r) = 42.4 Hz and uJs/C^Tr) = 120 Hz. The diamonds show the 
result found from using our approximation. 

Figure 6 The exact value for the specific heat is shown as the solid curve, and our 
approximation as the diamonds. = 2000 and the frequencies are as in the 
previous figure. 

Figure 7 The ratio of our approximate specific heat to the exact value is shown with 
diamonds. The result found from using the bulk result is shown with crosses. 

Figure 8 The ratio of the ground state occupation number found using approximation 
to the exact result is shown by the diamonds. The ratio of the bulk specific heat 
to the exact value is shown by the crosses. The frequencies and particle numbers 
are the same as the previous three figures. 
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